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Abstract 

A  non-isothermal  model  of  a  proton  exchange  membrane  (PEM)  fuel  cell  in  contact  with  interdigitated  gas  distributors  has  been  performed. 
The  model  accounts  for  the  major  transports  of  convective  and  diffusive  heat  and  mass  transfer,  electrode  kinetics,  and  potential  fields.  The 
effects  of  flow  orientation  and  total  overpotential  across  a  five-layer  membrane-electrode  assembly  on  the  thermal  behaviors  in  a  PEM  fuel 
cell  are  examined.  A  unique  feature  of  the  model  is  the  implementation  of  a  thermal-electrochemical  algorithm  to  predict  the  fluid-phase 
temperature  as  well  as  the  solid-matrix  temperature  in  a  PEM  fuel  cell.  The  simulation  results  reveal  both  the  solid-matrix  temperature  and  the 
fluid-phase  temperature  are  increased  with  increasing  total  overpotential.  Moreover,  the  fluid-phase  and  solid-matrix  temperature  distributions 
are  significantly  affected  by  the  flow  orientation  in  the  PEM  fuel  cell.  Replacing  the  parallel-flow  geometry  by  the  counter-flow  geometry  has 
an  advantage  of  reducing  the  local  maximum  temperature  inside  the  fuel  cell.  Thermal  effects  on  the  active  material  degradation  and  hence 
fuel  cell  durability  will  be  incorporated  in  the  future  work. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  irreversibility  of  electrochemical  reactions  is  the 
majority  of  heat  source  in  a  proton  exchange  membrane 
(PEM)  fuel  cell,  which  raises  the  fluid  temperature  as  well 
as  the  solid-matrix  temperature  during  operation.  The  local 
temperature  distribution  inside  a  fuel  cell  has  a  strong  impact 
on  the  fuel  cell  performance  since  it  affects  the  water-vapor 
distribution  by  means  of  condensation.  Insufficient  cooling 
may  result  in  local  hot  spots  and  thus  dehydrate,  shrink 
or  even  rupture  the  membrane.  In  addition,  the  kinetics  of 
electrochemical  reactions  directly  depends  on  temperature. 
Therefore,  as  far  as  the  reliability  and  durability  are  con¬ 
cerned,  an  efficient  thermal  management  of  a  PEM  fuel  cell 
becomes  crucial.  It  includes  not  only  the  removal  of  the  gen¬ 
erated  heat  from  inside  the  fuel  cell  to  the  outside  but  also 
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the  spatial  uniformity  of  temperature  distribution  that  avoids 
local  hot  spots. 

In  the  past  decade,  numerous  research  efforts  have  been 
devoted  to  develop  realistic  simulation  models.  Notable  early 
works  include  Bemardi  and  Verbrugge  [1,2]  and  Springer 
and  his  co-workers  [3,4].  These  models  are  one-dimensional 
and  only  account  for  diffusive  mass  transport  and  electro¬ 
chemical  kinetics.  Later  on,  several  two-dimensional  models 
were  developed  by  Nguyen  and  White  [5],  Fuller  and  New¬ 
man  [6],  and  Gurau  et  al.  [7].  Most  of  these  assume  some 
concentration  profile  of  reactant  species  along  the  channel 
except  for  [7],  which  accounts  directly  for  convective  mass 
transport.  Recently,  Shimpalee  and  Dutta  [8]  described  a 
steady  state,  isothermal,  three-dimensional,  and  single  phase 
PEM  fuel  cell  model.  It  used  a  commercial  code  to  solve 
the  complete  Navier-Stokes  equations.  Zhou  and  Liu  [9], 
Um  and  Wang  [10],  and  Hwang  et  al.  [11]  described  a  3- 
D  model  for  PEM  fuel  cell.  Their  results  agree  well  with 
their  own  experimental  observations.  It  is  interested  to  note 
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Nomenclature 


Ci 


Dt 

F 

hv 

l 

J 

kc 


kf 

ks 

M 

P 

Pi,  m 


Pi,  s 

R 

S 

T 

u 

x,y 


mole  concentration  of  the  species  i  (mol  m  3), 


p 

RT 


diffusivity  of  species  i  (m 2  s-1) 

Faraday’s  constant  (96487  C  mol-1) 
interfacial  heat  transfer  coefficient  (volumet¬ 


ric)  (Wm3  K-1) 


current  density  (Am-2) 

transfer  current  density  (Am-3) 

thermal  conductivity  of  the  solid  phase  in  the 

catalyst  layer  (W  K-1  m-1) 

thermal  conductivity  of  the  fluid  phase 

(W  K-1  m-1) 

thermal  conductivity  of  the  solid  phase 
(W  K-1  m-1) 

molecular  weight  (kg  mol-1) 
pressure  (Pa) 

possibility  of  the  electrolyte  in  the  connection 
of  the  catalyst  layer 

possibility  of  the  catalyst  in  the  connection  of 
the  catalyst  layer 

universal  gas  constant  (W  mol-1  K-1) 
source  terms  in  the  governing  equations 
temperature  (K) 
velocity  vectors  (ms-1) 
coordinate  system,  Fig.  1  (m) 


Greek  symbols 
a  symmetric  factor 

s  porosity  (gas  diffusion  layer) 

£c  porosity  of  the  catalyst  layer 

rjtot  total  overpotential  across  the  MEA  (V) 

k  permeability  (m2) 

fi  viscosity  (ms-2) 

p  density  (kg  m3) 

crm  ionic  conductivity  of  the  membrane  phase 

(£2-1  m1 ) 

<rs  electric  conductivity  of  the  catalyst  phase 
(^-1m-1) 
r  tortuosity 

0m  potential  of  the  ionic  conductor  (electrolyte 
phase)  (V) 

0s  potential  of  the  electric  conductor  (catalyst 
phase)  (V) 

vm  volume  fraction  of  the  ionic  conductor  (elec¬ 

trolyte  phase)  in  the  catalyst  layer 
us  volume  fraction  of  the  electric  conductor  (cat¬ 
alyst  phase)  in  the  catalyst  layer 
co  mass  fraction 


Subscript 
a  anode 

c  catalyst  or  cathode 


eff 

effective 

e 

energy 

f 

fluid 

HOR 

hydrogen  oxidation  reaction 

i 

species 

j 

electricity 

m 

membrane  phase,  momentum 

0 

exchange  current  density 

ORR 

oxygen  reduction  reaction 

ref 

reference 

s 

solid,  catalyst  phase 

T 

transfer  current 

tot 

total 

that  all  above  predicted  results  are  based  on  the  isothermal 
model. 

The  present  study  is  aimed  to  develop  a  non-isothermal 
model  to  examine  the  thermal-electrochemical  characteris¬ 
tics  inside  a  PEM  fuel  cell.  A  set  of  conservation  equations 
of  mass,  momentum,  energy,  species  and  charge  is  developed 
and  solved  numerically  with  proper  account  of  electrochem¬ 
ical  kinetics  and  fluid  dynamics.  The  Brinkman  extension 
to  Darcy  flow  describes  the  fluid  flow  characteristics  in  the 
porous  electrodes.  The  Stefan-Maxwell  correlations  along 
with  the  Bruggman  modification  illustrate  the  multi- species 
diffusion  in  the  porous  electrode.  A  two-equation  approach 
is  used  to  account  for  the  local  thermal  non-equilibrium 
between  the  solid  matrices  and  the  fluids  in  the  gas  dif¬ 
fusion  layers.  In  the  catalyst  layers,  the  heat  generation 
due  to  overpotential  heating  is  determined  from  the  macro¬ 
scopic  electrochemical  model.  Computational  fluid  dynamics 
(CFD)  methodology  is  employed  to  integrate  electrochemi¬ 
cal  processes  with  co-transports  of  multi-physics  in  the  PEM 
fuel  cell.  The  present  model  is  one  of  the  first  endeavors  to 
simultaneously  predict  the  solid-matrix  temperature  and  the 
fluid-phase  temperature  inside  a  PEM  fuel  cell.  It  enables  a 
comprehensive  understanding  of  the  mechanisms  responsi¬ 
ble  for  thermal  pathways.  Most  importantly,  the  possibility 
of  hot  spots  within  a  PEM  fuel  cell  has  been  successfully 
assessed.  The  results  obtained  by  the  present  model  would 
be  beneficial  for  the  design  of  thermal  management  of  a  low- 
temperature  fuel  cell.  It  is  also  helpful  in  the  further  accurate 
analyses  of  the  fuel-cell  thermal  performance  by  considering 
the  temperature-dependent  physical  properties  inside  a  fuel 
cell. 


2.  Model  description 

A  sectional  view  of  a  typical  PEM  fuel  cell  is  given  in 
Fig.  1.  The  model  domain  is  confined  within  a  five-layer 
MEA  (membrane-electrode  assembly),  i.e.,  two  gas  diffu¬ 
sion  layers  (GDLs),  two  catalyst  layers  (CLs),  and  a  proton 
exchange  membrane  (PEM).  Each  GDL  is  in  contact  with  an 
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Fig.  1.  Sectional  view  of  the  fuel  cell  module. 


interdigitated  gas  distributor,  which  has  an  inlet  channel,  a 
current  collector,  and  an  outlet  channel.  Humidified  hydro¬ 
gen  and  air  are  supplied  to  the  inlet  channels  of  the  anode  and 
cathode,  respectively.  In  the  anodic  catalyst  layer,  hydrogen  is 
consumed  to  form  protons  that  carry  the  ionic  current  to  the 
cathode.  In  the  cathodic  catalyst  layer,  the  electrochemical 
reaction  not  only  consumes  the  oxygen  but  also  produces  the 
water.  The  hydrogen  oxidation  reaction  (HOR)  in  the  anodic 
catalyst  layer  and  the  oxygen  reduction  reaction  (ORR)  in  the 
cathodic  catalyst  layer  are,  respectively,  expressed  as  by  the 


(a)  Counter  flow 


following  equations: 

±H2  -*  H+  +  e- 

0) 

\02  +  H+  +  e“  ±H20 

(2) 

Both  feeds  are  regarded  as  ideal  gases  and  are  transported 
through  diffusion  and  convection.  The  electrodes  are  treated 
as  homogeneous  porous  media  with  uniform  morphological 
properties  such  as  porosity  and  permeability.  All  gases  within 
each  of  the  electrodes  exist  as  a  continuous  phase.  In  addi- 


(b)  Parallel  flow 


Fig.  2.  Schematic  drawing  of  the  parallel-flow  geometry  and  the  counter-flow  geometry. 
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tion,  water  in  the  vapor  phase  is  considered  to  simplify  the 
model. 

In  the  present  study,  two  kinds  of  flow  orientations  are 
examined,  i.e.,  the  counter  flow  and  the  parallel  flow  (Fig.  2). 

2.7.  Governing  equations 

The  ionic  and  electronic  current  balances  in  the  PEM  and 
GDLs,  based  on  the  Ohm’s  law,  are,  respectively,  described 
as  the  following  equations: 

V(-<rmV0m)  =  0  (3) 

V(-asV0s)  =  0  (4) 

where  0m  and  0S  are  the  membrane-phase  potentials  and  the 
catalyst-phase  potentials,  respectively.  crm  and  crs  are  the  ionic 
and  electronic  conductivities  of  the  PEM  and  GDLs,  respec¬ 
tively. 

In  the  catalyst  layer  of  a  PEM  fuel  cell,  the  porous  matrix 
contains  two  kinds  of  solid  phases,  i.e.,  ionic  conductor  (elec¬ 
trolyte)  and  electronic  conductor  (catalyst).  A  potential  dif¬ 
ference  exists  between  the  catalyst  and  electrolyte  to  drive 
the  transfer  current  (jj),  keeping  the  electrochemical  reac¬ 
tion  continuously.  The  current  passes  through  catalyst  layer 


can  be  decomposed  two  parts,  i.e., 

i  =  is  +  im  (5) 

where  is  and  im  are  the  currents  flowing  through  the  cata¬ 
lyst  and  the  electrolyte,  respectively.  Since  the  electrodes  are 
electroneutral  everywhere,  there  is  no  charge-buildup  in  the 
catalyst  layers.  Thus,  the  charge  conservation  is 

V  •  i  =  0  (6) 

That  is 

V  •  is  =  -  V  •  im  (7) 

These  two  current  components  interact  through  electro¬ 
chemical  reactions.  The  electrons  are  transferred  to  the  cata¬ 
lyst  from  the  electrolyte  in  the  anodic  catalyst  layer,  and  vice 
versa  in  the  cathodic  catalyst  layer.  Application  of  Ohm’s  law 
to  equation  yields  the  current  conservation: 


^  ’  (  orSfeffV0s)  —  Sj 

(8) 

V  •  (  crm,effV0m)  —  Sj 

(9) 

where  the  sources  terms  Sj  and  —Sj  are  the  local  transfer  cur¬ 
rent  densities  corresponds  to  the  HOR  and  ORR  in  the  anode 
and  cathode,  respectively,  creating  and  consuming  protons. 


Table  1 


Porous-electrochemical  properties  of  the  present  modeled  fuel  cell 


Description 

Unit 

Value 

Exchange  current  density 

Anode,  j0,a 

Am-3 

5.0  X  103 

Cathode,  j0,c 

Am-3 

1.0  x  10“3 

Symmetric  factor 

Anode,  aa 

— 

1.0 

Cathode,  ac 

— 

0.5 

Porosity 

a-GDL,  c-GDL,  e 

— 

0.5 

a-CL,  c-CL,  sc 

— 

0.5 

Permeability 

a-GDL,  c-GDL,  k 

? 

m 

1.57  x  10“12 

a-GDL,  c-GDL,  kc 

2 

m 

1.57  x  10“12 

Tortuosity 

a-GDL,  c-GDL,  r 

— 

1.5 

a-GDL,  c-GDL,  rc 

— 

1.5 

Thermal  conductivity 

GDL,  kS  Q ff 

Wkg-1  K-1 

1.7 

PEM,  km 

Wkg-1  K"1 

0.5 

Anode  gas,  fcf>e ff 

Wkg-1  K"1 

0.182 

Cathode  gas,  fcf^ff 

Wkg-1  K"1 

0.051 

Electric  conductivity 

GDL  (electronic),  crs 

£2_1  m-1 

300 

PEM  (ionic),  am 

m-1 

14.4 

Inlet  pressure 

Anode,  /?in,a 

kPa 

1.013  x  105 

Cathode,  pin,a 

kPa 

1.013  x  105 

Mass  flow  rate 

Anode 

kgs-1 

1.19  x  10“6 

Cathode 

kgs-1 

2.98  x  10“5 

Species  mass  fraction  at  cathode  inlet  (saturated  air  at  STP) 

Oxygen,  coo2 

— 

21.7% 

Water,  wh2o,c 

— 

2.1% 

Nitrogen,  o;n9 

— 

77.2% 

Total 

— 

100% 

Species  mass  fraction  at  anode  inlet  (saturated  H2  at  STP) 

Hydrogen,  coun 

— 

76.5% 

Water,  ct»H2o,a 

— 

23.5% 

Total 

— 

100% 
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They  are  denoted  as  j’x.hor  and  j’x.orr?  respectively,  in  the 
following  discussion.  <xs,eff  and  crm,eff  are  the  effective  elec¬ 
tronic  and  ionic  conductivities  of  the  catalyst  and  electrolyte, 
respectively.  They  are  modeled  as 

^s,eff  =  ^s(l  -  Sc)  X  US  X  pijS  (10) 

^nxeff  —  £c)  x  x  Pi, m  (H) 


where  us  and  um  are  the  volume  fraction  of  the  catalyst  and 
electrolyte  in  the  catalyst  layer,  respectively,  pt ,s  and  p^m  are 
the  possibilities  of  the  catalyst  and  electrolyte  in  the  connec¬ 
tion  of  the  catalyst  layer,  respectively  [7,1 1].  It  is  noted  that 
only  a  long-range  connection  of  the  same  particles  stretch 
through  the  entire  catalyst  layer  ensures  good  conductivity. 

The  present  model  takes  into  account  two  species  in  the 
anode  (H2,  H2O),  and  three  in  the  cathode  (O2,  H2O,  N2). 
The  species  transports  based  on  the  Stefan-Maxwell  multi- 
component  diffusion  are  given  by  the  following  equations: 


(  " 

pn  •  =  V  •  i  poj,)  Dixn\j 

{  j=  1 


M 

Ml 


[S7coj  +  coj 


+  {Xj 


Pp 

-  COj) - 

p  . 


(12) 


The  effective  diffusivities  of  the  species  i  in  the  porous  elec¬ 
trode  follows  the  Bruggman  model  [12],  i.e., 

A\e  ff  =  eTDi  (13) 


The  source  terms  Si  represent  the  consumption  of  the  reac¬ 
tants  during  the  electrochemical  reaction.  It  becomes  Si  = 
jx,ORR^o 2/4F  and  Si  =  jx,HOR^H2/2T1  in  the  cathodic  and 
anodic  catalyst  layers,  respectively.  In  the  gas  diffusion  layer 
it  is  nothing.  According  to  the  Butler- Volmer  correlation 
[13],  the  relationship  among  the  local  transfer  current  density 
(jx),  the  reactant  concentrations  (cf),  and  phase  potentials  (0S, 
0m)  can  be  described  as  the  following  equation: 


•  m 

Jx,OOR  “  Jo,< 


c02 


c02,ref 


cu2o 

CH20,ref 


exp 


4  arF 


(0m  0s) 


exp 


RT 

4(1  —  ac)F 
RT 


(0m  0s) 


(14) 


•  V 

Jx,HOR  —  Jo, 


a 


cH2,ref 


exp 


4  a^F 


RT 


(0s  0m) 


(15) 


where  j0, c  and  j0?a  are  the  cathodic  and  anodic  exchange  cur¬ 
rent  densities,  respectively. 

The  fluid  flow  in  the  porous  media  is  described  by  the 
following  equations: 


(16) 

(17) 


where  p  is  the  density,  pi  the  viscosity,  u  the  velocity  vec¬ 
tor,  and  p  the  pressure.  The  source  term  in  the  momen¬ 
tum  equations  is  based  on  the  Darcy’s  law,  representing  an 
extra  drag  force  proportional  to  fluid  viscosity  and  velocity, 
and  inversely  proportional  to  the  permeability  of  a  porous 
medium,  i.e.,  Sm  =  — (/x//e) u,  where  k  is  the  permeability. 

As  for  the  energy  equations,  the  two-equation  model  is 
used  to  describe  the  thermal  behaviors  in  the  gas  diffusion 
layer.  The  energy  equations  for  fluid  and  solid  phases,  respec¬ 
tively,  are 

(pCp)fU  •  V7f  =  V  •  (kf,effV7f)  -  Se,G dl  (18) 

0  =  V  •  (ks,effVTs)  +  iSe,GDL  (19) 

The  source  terms  Se,GDL  =  —hv-(Ts  —  7f)  represent  the 
thermal  interaction  between  the  solid  matrices  and  the  flu¬ 
ids.  hv  is  the  interfacial  heat  transfer  coefficient  (volumetric) 
between  the  solid  matrices  and  the  reactants  fluid  in  porous 
medium  [14].  The  effective  thermal  conductivities  of  both 
phases  are,  respectively,  defined  as 

£s,eff  =  (1  -  o)h  (20) 

£f,eff  =  skf  (21) 


Fig.  3.  Flow  chart  of  the  numerical  simulation. 


pu  •  Vu  =  —V p  +  V  •  (pVu)  +  Sm 
V(pu)  =  0 
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(a)  Parallel  flow 


(b)  Counter  flow 
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Fig.  4.  Comparison  of  flow  velocity  vectors  between  the  parallel  and  counter  flows. 


where  ks  and  kf  are  the  thermal  conductivities  of  the  solid 
matrix  and  the  reactant  fluid,  respectively. 

In  the  catalyst  layer,  the  electrochemical  reaction  occurs 
at  the  interface  of  reactant  fluid  and  catalyst.  Physically,  the 
fluid  and  solid  phases  in  the  catalyst  layer  have  the  same 


temperatures,  i.e., 

(pCp)fU  •  V7f  =  V  •  (&c,effV7f)  +  Se,CL 


Tf  =  Ts 


(a)  Parallel  flow 


(b)  Counter  flow 
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Fig.  5.  Comparison  of  pressure  distribution  between  the  parallel  and  counter  flows. 
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In  the  present  model,  the  radiation  heat  flux  and  the  energy 
dissipation  due  to  Joule  heating  are  neglected.  Therefore,  the 
source  term  in  the  above  equation  can  be  represented  by  the 
overpotential  heating  by  the  electrochemical  activation,  i.e., 
Se,CL  =  jx,HOR  x  (0m  —  0s)in  the  anodic  catalyst  layer  and 
iSe,CL  =  J*t,orr  x  (0s  —  0m)  in  the  cathodic  catalyst  layer, 
respectively.  The  effective  thermal  conductivity  of  the  cat¬ 
alyst  layer  is  determined  by  the  following  equation  [12]: 


^c,eff  —  2k  c  + 


1 


s/(2kc  +  kf)  +  (1  -  e)/3 K 


(24) 


where  kc  is  the  weight-averaged  conductivity  between  the 
ionic  conductor  (such  as  Nation™)  and  the  electric  conductor 
(such  as  Pt/C). 


As  for  the  PEM,  a  typical  conduction  equation  is  employed 
to  describe  the  thermal  behavior  in  the  impermeable  material, 
i.e., 

V  •  (*mvrs)  =  0  (25) 


where  km  is  the  ionic  conductivity  of  the  PEM. 

2.2.  Boundary  conditions 

The  electrochemical  and  physical  properties  used  in  the 
calculation  are  given  in  Table  1 .  The  temperatures  for  both 
feeds  along  with  the  surfaces  of  the  gas  distributor  are 
fixed  at  298  K.  Both  outlets  of  the  module  have  an  ambi¬ 
ent  pressure.  The  anode  is  supplied  with  the  humidified 


Max:  314.159  Max:  314.159 


Min:  298  Min:  298 


Fig.  6.  Effect  of  total  overpotential  on  the  solid-phase  temperature  distributions  for  the  parallel-flow  geometry:  (a)  ij tot  =  0.45  V;  (b)  ijtot  =  0.65  V. 
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Max:360.939  Max:360.94 


Min:298  Min:298 


Fig.  6.  ( Continued ). 


hydrogen  of  mass  fractions  of  76.5/23.5%  for  H2/H2O.  The 
cathodic  side  feeds  with  the  saturated  air  of  21.7/2.1/77.2% 
for  O2/H2O/N2,  where  N2  is  considered  as  an  inert  gas  and 
serves  as  diluents.  Reactants  delivered  to  the  cathode  and 
anode  are  2.98  x  1CT5  and  1.19  x  10  6 kgs  ^respectively. 
The  difference  of  electronic-conductor  potential  (</>s)  between 
two  contact  surfaces  between  the  electrodes  and  gas  dis¬ 
tributors  represents  the  total  overpotential  (77 tot)  across  the 
five-layer  ME  A.  The  potential  at  the  contact  surfaces  between 
the  c-GDL  and  the  current  collector  is  arbitrarily  chosen 
to  be  zero,  while  the  total  overpotential  is  used  as  bound¬ 
ary  condition  at  the  anodic  current  collector.  For  the  rest 
of  the  boundaries  they  have  either  insulation  or  symmetry 
conditions. 


2.3.  Numerical  methods 

The  solutions  of  the  above  equations  are  obtained  with  a 
general  purpose  commercial  solver,  Femlab.  It  uses  the  Broy- 
den’s  method  with  an  FU-decomposition  pre-conditioner  to 
solve  the  non-linear  equations  iteratively.  To  reduce  conti¬ 
nuity  errors,  a  penalty  term  is  employed  for  pressure.  Thus, 
there  is  a  continuous  part  of  the  pressure  and  piecewise  con¬ 
stant  part  providing  and  extra  DOF  (degree  of  freedom)  for 
pressure  on  each  element.  It  uses  Newton-Raphson  iteration 
to  solve  the  close-coupled  groups  (velocity,  pressure,  temper¬ 
ature,  concentration  and  electricity)  and  uses  the  frontal  algo¬ 
rithm  (Gaussian  elimination)  to  solve  the  linearized  system 
of  equations  for  each  iteration.  In  the  computational  domain, 
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a  total  of  1828  nodes  and  8789  elements  were  used  (quadratic 
velocities  in  each  direction),  with  a  fine  mesh  throughout  the 
region.  The  program  gives  results  within  1%  of  each  other 
on  the  finest  meshes  used.  The  iterations  proceeded  until  the 
change  in  the  calculated  air  flow  rate  between  20  consecu¬ 
tive  iterations  was  less  than  0.1%.  In  order  to  achieve  that, 
the  necessary  number  of  iterations  varied  between  200  and 
350.  The  CPU  time  ranged  from  10  to  100  min  on  a  Pentium 
IV  PC  (2.8  GHz,  2  GB  RAM)  using  Windows  XP  operating 
system.  Fig.  3  shows  a  flow  chart  of  the  present  numerical 
modeling. 


3.  Results  and  discussion 

Fig.  4  compares  the  flow  velocity  distributions  between 
the  parallel-flow  geometry  and  the  counter-flow  geometry. 
The  point  marks  on  each  plot  mean  the  local  maximum  and 
minimum  velocities  in  the  module.  It  is  seen  that  the  fuel 
side  (anode)  has  a  lower  velocity  than  the  air  side  (cathode) 
since  it  has  a  smaller  stoichiometric  coefficient.  Results  also 
show  that  the  maximum  velocity  occurs  at  the  cathode  exit 
around  the  rib  surface.  The  flow  around  the  corners  formed  by 
the  symmetric  planes  (y  =  0,  and  2.0  mm)  and  the  interfaces 
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Fig.  7.  Solid-phase  temperature  distributions  of  the  counter-flow  geometry  for  ^ot  =  0.65  V. 
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between  the  PEM  and  catalyst  layers  (x  =  0.4  and  0.5  mm)  is 
nearly  stagnant.  Fig.  5  shows  a  comparison  of  the  pressure 
distribution  in  the  parallel-flow  geometry  and  counter-flow 
geometry.  Both  outlets  of  the  module  are  opened  to  the  ambi¬ 
ent,  which  serve  as  reference  pressures.  It  is  seen  that  the 
pressure  drops  across  the  cathode  is  significantly  higher  than 
that  across  the  anode,  typically  about  1 .5  and  0.15  kPa  for  the 
cathode  and  the  anode,  respectively.  This  is  because  a  higher 
fluid  density  is  accompanied  by  the  cathode  flow,  and  a  more 
fluid  should  be  driven  across  the  cathode  (Table  1).  In  gen¬ 


eral,  the  pressure  drop  characteristics  are  not  affected  by  the 
flow  direction  in  essential. 

3D  mappings  of  the  solid-phase  temperatures  for  the 
parallel-flow  geometry  and  the  counter-flow  geometry  are 
shown  in  Figs.  6  and  7,  respectively.  The  corresponding 
projects  showing  the  isothermal  contours  are  also  provided 
in  these  figures.  The  point  marks  shown  on  each  plot  indicate 
the  maximum  or  minimum  temperatures  in  the  module. 

Fig.  6(a)  and  (b)  shows  the  effect  of  total  overpotential 
(rj tot)  on  the  solid-phase  temperature  distributions  for  the 


Fig.  8.  Comparison  of  fluid-phase  temperature  distributions  at  several  y  stations  between  parallel  flows  and  counter  flows,  rjtot  =  0.65  V. 
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Fig.  9.  Effect  of  total  overpotential  on  the  concentration  distribution  in  the  computational  module. 


parallel-flow  geometry.  It  is  seen  that  the  solid-phase  tem¬ 
perature  in  the  cathode  is  higher  than  that  in  the  anode. 
This  is  because  the  c-CL  dissipates  more  heat  by  the  electro¬ 
chemical  reaction  than  the  a-CL.  The  maximum  temperature 
occurs  in  the  c-CL  cutting  across  the  upper  symmetric  plane 
(y  =  2.0  mm).  In  addition,  the  region  near  the  rib  surface  has 
a  low  solid  phase  temperature.  Near  the  fuel  and  oxidant 
entrances,  the  solid  phase  has  low  temperatures  due  to  the 
significant  forced  convection  by  the  inlet  fluid.  When  the  total 
overpotential  increases  form  ijtot  =  0.45-0.65  V,  as  shown  in 
Fig.  6(b),  the  solid-phase  temperature  distribution  becomes 
more  uneven.  In  addition,  the  maximum  solid-phase  temper¬ 
ature  increases  from  314  to  361  K. 

Fig.  7  shows  the  results  of  the  counter-flow  geometry  with 
the  same  total  overpotential  as  that  in  Fig.  6(b).  Again,  the 
minimum  solid-phase  temperatures  occur  in  the  region  near 
the  rib  surfaces  (y  =  0  and  0.9  mm).  However,  the  maximum 
solid-phase  temperature  has  moved  to  the  bottom  of  the  mod¬ 
ule,  i.e.,  at  the  c-CL  cutting  across  the  middle  of  cathode  exit 
(y  =  0).  The  peak  of  the  solid-phase  temperature  is  reduced 
from  361  to  343  K  with  the  parallel-flow  geometry  instead 
of  the  counter-flow  geometry.  This  is  because  the  heat  dissi¬ 
pated  by  the  overpotential  heating  near  the  cathode  outlet  is 
somewhat  cooled  down  by  the  anode  inlet  flow  that  removes 
the  conduction  heat  through  the  PEM. 

Fig.  8(a)  and  (b)  shows  the  fluid-phase  temperature  dis¬ 
tribution  along  several  elevations  (i.e.,  y  =  0.2,  0.6,  1.0,  1.4, 
1.8  mm)  cutting  across  the  module  for  the  parallel-flow  geom¬ 
etry  and  the  counter-flow  geometry,  respectively.  Since  the 
PEM  is  impermeable,  no  data  are  shown  in  the  region  of 
0.4  mm  <  v  <  0.5  mm.  In  both  flow  geometries,  the  fluid-phase 
temperature  in  the  cathode  is  higher  than  that  in  the  anode. 
This  is  because  the  significant  heat  generation  in  the  cathodic 


catalyst  layer.  A  stronger  electrochemical  kinetics  of  HOR 
(i.e.,  higher  exchange  current  density,  Table  1)  in  the  a-CL 
requires  a  less  overpotential  to  drive  the  through-flow  current 
in  the  fuel  cell.  Consequently,  a  less  overpotential  heating  in 
the  a-CL  results  in  a  lower  temperature.  The  fluid-phase  tem¬ 
perature  increases  downstream  due  to  the  heat  accumulation 
for  both  electrodes.  The  maximum  fluid-phase  temperature 
for  the  parallel-flow  geometry  is  higher  than  the  counter-flow 
geometry. 

Fig.  9  shows  the  effect  of  total  overpotential  on  the  con¬ 
centration  distribution  in  the  module  for  the  parallel-flow 
geometry.  The  data  shown  in  the  plots  are  normalized  by 


x/m 

Fig.  10.  Voltage  of  the  catalyst  phase  and  membrane  phase  along  the  mid¬ 
plane  of  the  module. 
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the  corresponding  inlet  concentration,  i.e.,  cou2  /  con2  for 
the  anode  and  coo2/coo2jm  for  the  cathode.  It  is  seen  from 
these  figures  that  an  increase  of  total  overpotential  increases 
the  fuel  and  oxidant  utilization  in  both  electrodes.  At  a 
fixed  total  overpotential,  the  local  concentration  decreases 
along  the  flow  direction.  Also,  it  decreases  as  the  flow 
approaches  the  catalyst  layer.  It  is  further  seen  that  the  con¬ 
centration  gradient  along  the  v  direction  in  the  cathode  is 
more  significant  than  that  in  the  anode.  This  is  because 
cathode  has  a  higher  velocity  (Fig.  4)  that  dominates  the 
species  transport  via  forced  convection.  In  contrast,  the  flow 
velocity  in  the  anode  side  is  relatively  stagnant  (Fig.  4). 
The  diffusion  becomes  more  significant  in  the  species 
transports. 

Fig.  10  shows  the  distribution  of  the  phase  potential  at  the 
elevation  cutting  across  the  module  mid-plane  (y=  1.0  mm). 
The  solid  lines  are  the  catalyst-phase  potential  while  the 
dashed  line  is  the  membrane-phase  potential.  The  differences 
between  the  above  potentials  in  the  catalyst  layers  represent 
the  activation  overpotentials,  i.e.,  =  </>c  —  0m  in  the  a-CL 


and  r)c  =  (f>m  —  0C  in  the  c-CL.  The  activation  overpotentials 
in  both  catalyst  layers  increases  along  the  depth  of  the  cata¬ 
lyst  layer.  That  is  the  largest  activation  overpotentials  occur 
at  the  both  sides  of  the  PEM.  The  activation  overpotential 
drop  in  the  c-CL  is  significantly  higher  than  that  in  the  a-CL. 
In  the  GDLs  and  PEM,  a  linear  drop  of  phase  potential  is 
found.  The  slope  in  the  PEM  is  higher  than  that  in  the  GDLs 
meaning  that  the  Ohmic  loss  in  the  PEM  is  higher  that  in  the 
GDLs. 

Fig.  11(a)  shows  a  full-field  distribution  of  the  local  cur¬ 
rent  density  in  the  module.  It  is  seen  that  the  current  directs 
form  the  surfaces  in  contact  with  the  anodic  gas  distributors 
toward  the  surfaces  in  contact  with  the  cathodic  gas  distribu¬ 
tors.  The  current  density  is  high  near  the  rib  surface  corner  and 
is  low  near  the  module  corner.  Fig.  11(b)  is  a  local  magnifi¬ 
cation  showing  the  current  density  distribution  in  the  catalyst 
layers  near  the  module  middle.  As  shown  in  Fig.  11(b),  the 
current  density  in  the  a-CL  gradually  decreases  along  the 
current  direction,  i.e.,  from  the  interface  of  the  GDL  and  the 
catalyst  layer  (i.e.,  v  =  0.35  mm)  to  the  PEM.  This  is  because 


x  io-3 


Fig.  11.  Distributions  of  current  density  vectors:  (a)  the  entire  computational  module;  (b)  magnification  around  the  catalyst  layers. 
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of  the  HOR  consumes  the  current.  In  contrast,  the  current 
density  is  gradually  recovered  in  the  c-CL  from  the  interface 
of  the  PEM  and  the  c-CL  (x  =  0.50  mm)  to  the  c-GDL  due  to 
the  ORR. 


4.  Concluding  remarks 

A  multi-physics  model  has  been  developed  to  investigate 
the  thermal-electrochemical  transports  in  a  PEM  fuel 
cell.  Conservative  equations  governing  the  co-transports 
of  mass/momentum/heat/species/charge  are  numerically 
solved  with  proper  account  of  electrochemical  kinetics.  In 
the  catalyst  layer,  a  general  energy  equation  is  derived  using 
the  volume-averaging  technique,  along  with  a  local  heat 
generation  resulting  from  electrochemical  reactions.  In  the 
gas  diffusion  layer,  a  two-equation  thermal  transport  model 
is  developed  to  resolve  the  fluid  and  solid  phase  temperatures 
for  the  first  time.  Results  show  that  both  the  solid-matrix 
temperature  and  the  fluid-phase  temperature  increase  with 
increasing  the  total  overpotential  (%*).  Under  the  same 
total  overpotential,  the  maximum  solid-phase  temperature 
is  reduced  by  replacing  the  parallel-flow  geometry  with  the 
counter-flow  geometry. 

The  present  paper  has  provided  an  innovative  aspect  in 
the  heat  transfer  of  fuel-cell  related  studies.  It  has  success¬ 
fully  predicted  the  fluid  and  solid  phase  temperatures  inside 
a  fuel  cell  simultaneously  for  the  first  time.  It  is  important 
to  understand  the  electrochemical/thermal  coupled  mecha¬ 
nisms  responsible  for  thermal  pathways  in  a  low-temperature 
fuel  cell.  The  thermal  effect  on  the  active  material  degra¬ 
dation  and  hence  fuel  cell  durability  will  be  studied  in  the 
future  work. 
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